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Abstract 

We present a novel approach for model reduction of nonlinear dynamical systems based on proper 
orthogonal decomposition (POD). Our method, derived from Density Matrix Renormalization 
Group (DMRG), provides a significant reduction in computational effort for the calculation of 
the reduced system, compared to a POD. The efficiency of the algorithm is tested on the one 
dimensional Burgers equations and a one dimensional equation of the Fisher type as nonlinear 
model systems. 
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Nonlinear dynamical systems arise in many fields as physics, e.g. turbulence math- 
ematics and biology |2|. They often require a significant high number of degrees of free- 
dom(dof) for simulation with suitable accuracy. For a large class of systems the solutions 
are regular, nevertheless the influence of the nonlinearity is essential. Here model reduction 
(MR) can lead to an efficient description, if the dynamics is effectively confined to a lower 
dimensional attractor in phase space. 

The aim of this work is to develop an algorithm that can find a reduced model for a given 
system, usually derived from a partial differential equation. 

One method to obtain such a reduced description is the so called proper orthogonal 
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decomposition |3| (POD). It is obtained by calculating the eigenvectors of the spatial 
covariance matrix of the field over the phase space, typically by using a empirical spatial 
covariance matrix from a number of realization of the dynamic evolution. The method itself 
is linear in that the phase space is reduced to a subspace (in which the relevant Attractor 
has to be embedded). Nevertheless it accounts for the nonlinearity and gives the 'optimal' 
linear reduction possible. By definition the POD requires a simulation of the full, unreduced 
system and a diagonalization of a symmetric matrix of similar size. While the later can be 
circumvented by the method of snapshots of Sirovich [Ll, the simulation in unavoidable. 

Within our approach, we try to calculate 'approximate' POD modes without simulating 
the full, unreduced system. This is achieved by following concepts from density matrix 
renormalization. Effectively, we adapt a system of already low dimensionality, to reproduce 
the full dynamics optimally. Consequently, all calculations are performed on low dimensional 
systems. In exchange several calculation steps have to be performed, but their number is 
proportional to the size of the full system of interest. Practically this is one possible way to 
study large dynamical systems, although within this work we are still restricted to spatially 
one dimensional systems. Beside from possible benefits for the efficiency, an interesting 
question is whether it is possible to reconstruct dynamic behavior of a system from the 
knowledge of subsystems only. 

This paper is organized as follows. First, we introduce the formulation of the equations 
defining the dynamical system. This includes the use of higher order tensors to describe the 
nonlinear part of the generator of the time evolution. Further, the discretization of the three 
model equations, namely the linear diffusion equation, the Burgers equation and a nonlinear 
diffusion equation, is presented. Then the type of orthogonal projection, which is used in 
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this work, is introduced and the basic concept of the proper orthogonal decomposition is 
recapitulated. After a brief outline of the single particle DMRG approach, the method 
devised in this paper is presented. The numerical results, including a comparison of our 
method with standard techniques, are given consecutively. This is followed by a short 
discussion of our approach and the conclusions. The appendix finally contains an analysis 
of the optimal reduction for the linear case. 

I. THE PROBLEM 

A. The Dynamical System 

We consider discretized versions of nonlinear evolution equations of the form 



Here $ is the field, G'($) is the nonlinear generator of evolution and we make use of the 
sum convention. The contributions L, Q and K are the linear, the quadratic and the cubic 
part, respectively, of the generator of evolution. Higher order terms can also be considered, 
but the number of nonlinear terms should be finite for our approach. Note that Q and K 
are third and fourth order tensors and have the corresponding transformation properties. 
The dynamical system described in Eq.([T]) is typically derived from a partial differential 
equation (PDE). The spatial discretization is then done by finite differences or equivalently 
by linear finite elements(FE). The temporal discretization is done by the simple explicit Euler 
method, although this choice is not relevant for our method. Here, we restrict ourselves to 
the spatially one dimensional case. In the following we exemplify our approach on simple 
toy problems. 

B. The Linear Diffusion Equation 

The diffusion equation describes diffusive transport of a scalar field, e.g. heat transport, 
in a medium. For homogenous media it is given by 
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with the diffusion constant d. Spatial discretization of the interval [0, 1] with nodes gives 
for the discrete Laplace operator with second order accuracy in Ax the following N x N 
matrix 
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Here homogenous Neumann conditions are assumed for a; = and x = 1, the spatial 
discretization step size is Ax = The explicit Euler method gives for the discrete time 
evolution with time step size ht the following equation 

<I(X, tn+l) = tn) + dhtAN^{x~tn) (4) 

where $ and x are iV-dimensional vectors, indicated by ~. Thus the linear part L in Eq. ([1]) 
is given by 

L = dhtAN. (5) 
The nonlinear contributions in Eq.([T]) vanish for the linear diffusion equation. 



C. The Burgers Equation 



As one nonlinear example we consider the Burgers equation [sj. It 
as well as a convective transport of a scalar field ^ and is given by 

d 

— $ = dA^ + z/($V)$. 



describes a diffusive 
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This equation is similar to the linear diffusion equation Eq.([2]) but with an additional term 
z/($V)$, describing the convection. This term is quadratic in the field $ and can be dis- 
cretized in the form of Q in Eq.([T]). For one space dimension, the V operator is simply the 
spatial derivative. This can be discretized with second order accuracy in Ax as 
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The term ($V) is also known as convective derivative. In ID the discretization is given by 
multiplying the rows of D^^n with the components of $: 

($V)iv,*j = ^iD^^N,i,j (8) 
here z, j indicate the component of the matrix/vector. Choosing 

Qi,j,k '■= ^Dx,N,j,kSij (9) 

gives a discretization of the convection term, as defined in Eq.(l8]) 

j,k j,k 

= iyYl ^iD.,N,^,k<^k = y ($V)^ <!>. (10) 

k 

D. Nonlinear Diffusion 



We consider here a Diffusion equation with a nonlinearity_^tliat resembles the action- 
equation. In particular 



potential part of the one dimensional FitzHugh-Nagumo(FN) [3, ^ 
the dynamics is defined by 

d 

= A<I)-<I)(l-<l>)(a-$) (11) 

where a is a constant. Eq.l ITT]) has stable equilibria at $ = and $ = 1 and an instable 
equilibrium at $ = a. The nonlinear term is cubic in the field. It can be rewritten as 
— $(1 — $)(a — $) = — + (1 + a)$^ — a$. Here the powers of $ are defined component 
wise. The cubic part — e.g. is discretized by 

= -SijSikSii (12) 

since 

J2 i-^iAk^ii) ^I'^k^i = (13) 

j,k,l 

Similarly, the quadratic part becomes 

Qi,j,k = (1 + a) 6ij6ik (14) 

and the linear part together with the contribution from the diffusive term is 

Li J = dhtAN,i,j - a6ij. (15) 

Since we deal only with discretized fields $ in the following, we drop the notation ~ for 
discrete variables. 
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E. The Reduction 



To obtain a reduced model, we will project the phase space to a lower dimensional 
subspace. Thus the reduction is linear, which simplifies the problem significantly. For 
nonlinear reduction see e.g. 

relevant subspace, e.g given by the column vectors of a x M matrix B, where is the 
dimensionality of the phase space and M that of the subspace (A^ > M). We will always 
assume an orthonormal basis in the following since B can always be brought to this form. 
This basis spans the range of the projection operator P, which is defined by 

P = BBK (16) 

The reduced dynamics is given by 

dtP^ = PLP^ + P (QP$P$) + P (KP$P$P$) (17) 

One can write this equation directly for the reduced phase space which is only M dimensional 
by using 

I) = P"f$ (18) 
L = B^LB (19) 

Qi,j,k = B\^^Qa,b,cBb,jBc^k (20) 
a,b,c 

^i,j,k,l = ^ Bl j<a,b,c,dBbjBc^kBd,l (21) 
a,b,c,d 

This gives the reduced equations which are still in the form of Eq. ([1]) as 

dt^ = Lij^j + Qijk^j^k + K^jkl^j^k^l (22) 

While this dynamics is defined on a smaller, M dimensional phase space it has to be noted 
that the operators are now typically dense, i.e. most entries in the tensors L, Q and K 
are nonzero. To define the reduction we have to make a choice for the relevant degrees of 
freedom that span the range of P, i.e. the orthonormal basis (ONB) B. Natural criteria for 
the determination of B would be based on the difference 

E(t) = <l>(t) - BB^^{t) = (1 - P)$(t), (23) 
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FIG. 1: Assembly of the superblock Laplace operator 
Initialisation: 
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FIG. 2: Graphical illustration of the DMRG initialization (or warmup) scheme 



e.g. the -error 

i?(t):=||E(t)||2. (24) 

For linear systems, as e.g. the linear diffusion equation, it can be shown (see section IVl|l 
that the projector onto the eigenstates with lowest absolute eigenvalue of the generator of 
the time evolution (i.e. L in Eq.([T])) lead to a minimal L^-error for long and short enough 
times. In the long time limit this approximation becomes even arbitrary accurate as long as 
the discarded eigenvalues are smaller than zero. This can be extended to nonlinear systems 
using the Proper Orthogonal Decomposition(POD). 
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FIG. 3: Graphical illustration of the DMRG iteration (or sweeping) scheme 



II. PROPER ORTHOGONAL DECOMPOSITION 



The proper orthogonal decomposition is a linear projection method which is widely usee 
in model reduction. On this topic an extensive literature exist. Some examples are 

m 

A short explanation of POD together with the method of snapshots is also given 
in 15|. One of the advantages of this method is the possibility to incorporate information 
from the nonlinear dynamics to obtain a linear reduction. The basic idea is to generate 
sample trajectories by simulating the dynamical system of interest. Then the spatial two 
point correlation matrix C is calculated 

:= (25) 

Here ()rp denotes an average over all sample trajectories. The eigenvectors of this symmetric 
matrix, which correspond to the highest eigenvalues span an 'optimal' subspace in the sense 
that the average least square truncation error 

e:={\mx,t)-P^x,t)\W (26) 

is minimal, see e.g. . The (orthonormal) basis vectors of this subspace constitute the 

columns of the matrix B which defines the projection operator P, see Eq. (fT6l) . The prac- 
tical application in the following algorithm is simple: Once the operators L, Q, and K are 
calculated and the initial conditions are given, we can simulate the dynamics of the field 
e.g. by Eq.([T]). For the reduced system Eq. (l22|) is used instead. During the simulation the 
data for the covariance matrix C is accumulated, if necessary several simulation runs are 
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performed using different initial condition with appropriate weigliting. The eigenvectors of 

n 

C are calculated using standard methods [6|. Then B is constructed from those eigenvec- 
tors corresponding to the highest eigenvalues. The discarded eigenvalues can also provide 
information on the quality of the reduction. The POD can be applied relatively independent 
from the actual blocking method. It should be noted that such a reduction is optimal for 
describing the generated dataset, but not necessarily optimal in reproducing the underlying 



dynamics 
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III. BLOCKING METHOD 



Blocking methods were considered already earlier, e.g. [18|, mainly because in many 
problems not all spatial regions are of similar interest. Our motivation is different, we aim to 
decompose a calculation into more feasible sub-problems. In the linear case the basis B can 
be calculated in principle by simply diagonalizing L. Technically this is the same problem as 
determining the eigenstates of the Laplace operator, which also describes the single quantum 



mechanical ID-particle in a box. S.White 
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20| used this toy model for introducing the 



so called Density Matrix Renormalization Group (DMRG). This approach has been since 
then applied most successfully to quantum many-body problems. In the following we want 
to carry this analogy further. Instead of approximating the eigenvalues/states of a linear 
operator we use a similar algorithm to obtain an approximate POD of a nonlinear system. 
To this end a few modifications are necessary, so we first summarize the original DMRG 
method. 



A. Single Particle DMRG 

In the DMRG toy problem the system is split into blocks of size m. For each block a 
block-Laplace operator is stored, as well as the links T that define the interaction with the 
neighboring sites. The assembly of the superblock Laplace operator is pictorially presented 
in Fig.lll This superblock operator is a (2m+2) x (2m+2) matrix and has to be diagonalized. 
From the eigenstates the so called target states 0*, i = 1 . . . (m + 1) are selected, usually the 
low lying spectrum. Since we have in this special case a single particle problem, the block 
truncation matrix R can be calculated simply by applying a Gram-Schmidt orthogonalization 
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to the block part of the target states. 

Rij = Gram-Schmidt (0}) i = 1 . . . {M + 1), j E Block (27) 

For more general DMRG applications R would be calculated by diagonalizing the density 
matrix of the target state, i? is a (m + 1) x m matrix which projects the phase space of one 
system side to an effective block. 

The effective block Laplacian Le// and the effective block link Tejf are derived as follows 

^eff = R^ LsideR Teff = R^Tside (28) 

where the form of Lgide and Tgide are depicted in Fig. [1] for the left side. By this process 
effective blocks are calculated, that describe a higher number of sites, but have numerically 
still m degrees of freedom. 

In DMRG this 'growing' of blocks is first used in an initialization step until the superblock 
describes a sufficiently large system, see Fig. [21 Then an iteration is carried out to increase 
accuracy. Here only one side is grown, while the other is replaced by an already calculated 
block so that the effective size of the superblock is constant, see Fig|3l 

B. DMRG-POD 

Basically three modifications are necessary to obtain a DMRG-POD algorithm 

First, instead of a diagonalization of the superblock operator, a POD on the superblock 
system has to be performed. This is composed of first, a simulation of the superblock 
system, as defined in Eq.( !22l) . Then the superblock correlation matrix from the gener- 
ated data has to be diagonalized. This gives an orthonormal set of vectors which are 
the target vectors in the context of DMRG. 

Second, to each sub-block there exist not only a linear sub-block operator but also higher 
order operators, given by third and higher order tensors, see Eg. (11 71) . These have to 
be updated in a similar way. 

Third, for the POD the initial states for the sample trajectories are crucial. The initial 
states are defined for the full system. They have to be projected onto the superblock 
system which requires all truncation matrices explicitly. 
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Concerning the first point, this is no great difference, since the POD (simulation and diago- 
nahzation of the correlation matrix) returns also a orthonormal set of 'relevant' states (the 
POD modes) that serve as target states, as described above. 

Beside the linear operator (L in Eq.([T])) which is assembled identically as the superblock 
operator in DMRG, the higher order operators have to be a assembled also. This is princi- 
pally possible, but complex. Here we use a simple trick. For all our models it is sufficient to 
know the component-wise squaring operator := ^ijSik- (And in some cases the deriva- 
tive operator which is linear and is also assembled like the superblock Laplace operator.) 
Q is purely diagonal, so no links have to be stored and assembled. The reduction with a 
truncation matrix R is straightforward: 

^i,j,k ~ ^ ^ Ri,a^a,b,cRb,jRc,k (2"^) 
a,b,c 

From this the higher order tensors can be calculated directly, e.g. for the Burgers equation 

^Burgers,ij,fc ■= T^^^ij^ijDx,N,i,k 
I 

= ''^^^j,i,lDx,N,l,k- (30) 

I 

For fourth and higher order operators this procedure is a bit memory consuming. E.g for 
calculating $^ it is more efficient to calculate first $tmp := = and then <l>^ = <l^<l = 

The third point is a small disadvantage, since the projection operators have all to be 
stored, which is not necessary in DMRG if only the energy values are of interest. However, 
here as well as in DMRG it is possible to expand a superblock state to a state of the original 
system as well as project down a system state to the superblock if all truncation matrices are 
stored. The down-projection of the A^- dimensional state is in particular done by iteratively 
contracting the m + 1 outermost sites of e.g. $ with the corresponding block truncation 
matrix R. Apart from the memory requirement this is simply a book keeping problem. 

It should be noted that only m + 1 most relevant states from the POD are used as target 
states. Thus only m+1 relevant states of the superblock are optimized although it represents 
2m + 2 dofs. This has to be considered in comparing the results. However, the DMRG POD 
is nevertheless faster than the full POD, see section IV A[ 

To summarize: Apart from the POD itself, which is a standard technique, no fundamental 
changes have to be implemented to get a DMRG-POD method from the simple toy model 
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DMRG. The assembly of linear operators has to be performed in any case, only our method 
requires several operators. The assembly of the Q operator is even more simple, since all links 
vanish. Reconstruction of full system states is also possible in DMRG, only it is mandatory 
for our method for evaluating the correct initial conditions. 

IV. APPLICATION 

For all applications we choose a finite differencing scheme of second order accuracy, ho- 
mogenous Neumann conditions at the boundaries and the explicit Euler method for the 
calculations. The details are given in section [H above. The boundary conditions as well 
as the time integration method can be chosen - more or less - arbitrarily. However, higher 
order finite elements in the spatial discretization lead to additional interactions between 
single dofs, i.e. a form of non-locality and do thereby complicate the problem. For the 
reduced system size always four dofs were retained. This is mainly for convenience and easy 
comparison. The success of the method does not depend strongly on this choice. 

As explained above, we measure the quality of a reduction by the L^-error, see Eq. (!24l) . 
It has the same units as the fields $ which are not further specified. The time units are also 
arbitrary. The error calculations in the following are performed in a separate program which 
gets the optimized bases from the various methods as input. Thus the simulation time do 
not have to coincide with the length of the POD simulations. Further we have chosen a 
different random seed for statistical initial conditions unless otherwise stated. 

A. Linear Diffusion 

For this problem the dynamics is given by Eq.([2]). The only nonzero contribution 
according to Eq.([T]) is L = A/^r. The eigenstates of L are the sine/cosine or Fourier modes 
whose contributions decay over time with characteristic life-time inversely proportional to 
the frequency /energy. Standard DMRG can be viewed as an approximate diagonalization 
method for an linear operator. Therefore it is very effective to find the optimal reduction 
determined by the eigenstates, see Appendix IVIl in the linear case. In contrast to the 
diagonalization, POD as well as our method depends on the initial conditions for the sample 
trajectories over which the averaging is carried out. Both POD approaches cannot exploit 
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FIG. 4: Reduced diffusion equation L^-Error E{t) for the analytical reduction (Fourier modes) the 
full POD and DMRG POD after initialization and several iterations, statistical initial condition, 
N=40, hf = 0.001 and d = 0.05. The error is expressed in units of ^, for the time axis arbitrary 
units are employed. Note that for clarity not all data points are shown as symbols. 
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FIG. 5: Reduced diffusion equation L^-Error E(t), identical statistical initialization. The error is 
expressed in units of <I>, for the time axis arbitrary units are employed. Note that for clarity not 
all data points are shown as symbols. 



the linearity of the evolution equation. This affects the quality of the results for linear 
problems compared to diagonalization-based methods. Nevertheless, restriction to a few 
sample trajectories can also be an advantage, since sometimes the interest lies on a certain 
region in phase space. However, for the diffusion equation we choose normally distributed 
initial conditions, i.e. the field ^o{xi) is normally distributed. This is then also true for 
the Fourier modes. By this choice effectively the whole phase space will be sampled for 
a high enough number of realizations. This is also due to the invariance of Eq.([2]) under 
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multiplication with a constant factor. 

For the POD it is important to integrate over long enough times. For short times the 
state moves in the direction of the highest frequency modes which are decaying most rapidly. 
Thus POD would give the wrong relevant modes. The POD is in fact not a very appropriate 
tool to reduce the whole phase space of the diffusion equation. In Fig. H] the error of the 
reduced fields $ is plotted in dependence of time. There the time step was dt = 10~^ and the 
diffusion constant d = 0.05. The spatial resolution was 40 lattice sites within the interval 
[0, 1]. In each POD step as well as for the error calculation the ensemble average, compare 
with Eq.( l25i) . has been averaged over 50 realizations of the initial conditions. From this 
result we can state several things. First, all POD-based methods show a remaining error in 
the long time limit. Second, the initialization steps of DMRG POD gives already reasonable 
results. An improvement due to the iteration is present, too. Third, our algorithm is able to 
compute the optimal reduction with even higher accuracy than the full POD. The last point 
is only paradox on the first glance. The inaccuracy of the full POD is in this case influenced 
from the statistical initial conditions, in order to sample the full phase space. Within our 
algorithm, much more initial conditions are taken into account as the superblock POD is 
performed repeatedly. This leads to a better statistics. In Fig. [5] we have shown the same 
results but using always the same initialization for calculating all PODs (but of course not 
for the error calculation). It is clear that in this case, our method has no advantage over the 
full POD anymore. On the other hand, the results from our algorithm are not worse than 
that from the full system POD, which is not clear a priori. 

B. Burgers Equation 

The Burgers equation is given by Eq.(l6]). The discretization used here is already described 
above. We begin our analysis with the choice of deterministic initial condition for the 
calculation of all PODs. In particular it is of the form 

$(f = 0, Xi) = e-^°("»-i)' Xi = 0...1. (31) 

Fig. [6] and [7| show the results for the L^-Error of the evolution. Here we have used two spatial 
resolutions, i.e. = 40 and = 100 nodes. The results are very similar. In contrast to the 
previous calculations the simulation runs for the error calculation are longer than the POD 
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FIG. 6: Reduced Burgers equation L^-Error E{t), deterministic initial condition, N=40, ht = 0.02, 
d = 0.01 and u = 0.1. The error is expressed in units of <I>, for the time axis arbitrary units are 
employed. Note that for clarity not all data points are shown as symbols. 
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FIG. 7: Reduced Burgers equation L^-Error E{t), deterministic initial condition, N=100, ht = 
0.005, d = 0.01 and z/ = 0.1. The inset shows the begin of the error evolution enlarged. The error 
is expressed in units of for the time axis arbitrary units are employed. Note that for clarity not 
all data points are shown as symbols. 

runs. The vertical line indicates the time interval of the POD runs. Here we have to state 
that the Fourier mode reduction is not optimal, which is not surprising as we consider a 
nonlinear system and a very particular region of phase space. Further, we see that the error 
curves show a very pronounced minimum after which the approximation seemingly breaks 
down. The corresponding time point lies well after the POD time-span. These minima 
correspond to the fact that after the passing of the wavefront the profile becomes fiat. The 
approximations do not reproduce the average value accurately, but show a spurious drift. 
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FIG. 8: Reduced Burgers equation L^-Error E{t), statistical initial condition, N=20, ht = 0.01, 
d = 0.05 and u = 0.1. The error is expressed in units of <I>, for the time axis arbitrary units are 
employed. Note that for clarity not all data points are shown as symbols. 

The passing of the reduced (flat) states by the original (flat) state creates the minima in 



It is remarkable that our methods yield better results than the POD within the POD 
time, even for the initialization. Here it should be recalled that the POD is optimal only 
for reconstructing the states used in the calculation. As stated above, the reconstruction of 
the dynamics that created these states, is a different thing as can be directly seen from our 
results. 

We continue our analysis of the Burgers equation by considering statistical initial condi- 
tions. In contrast to the calculations for the diffusion equation we have only three randomly 
sampled parameters in the initial condition. It is given by a peak of various height H, width 
W and position X. In particular it is deflned by the following equation 



Here, H and W are normally distributed whereas X is uniformly distributed. 

The results are shown in Fig. [HI For all methods the error reaches a plateau very quickly. 
The performance of the full system POD is slightly better than that of the DMRG POD. 
However, the errors from our approach are of the same order as from the full POD and one 
magnitude better than that of the Fourier mode based reduction. Also the iteration brings 
an improvement which reaches saturation already after the flrst step. 

For deterministic initial conditions the evolution of the error is not monotonic in contrast 
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FIG. 9: Reduced nonlinear diffusion equation L^-Error E{t), statistical initial condition, N=30, 
ht = 0.03, d = 0.01 and a = 0.5. The error is expressed in units of ^, for the time axis arbitrary 
units are employed. Note that for clarity not all data points are shown as symbols. 
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FIG. 10: Computing time for various system sizes and approaches, obtained by the reduced Burgers 
equation, statistical initial condition, N=40, ht = 0.005, d = 0.01 and ly = 0.1. 



to the case of statistical initial conditions. This is due to the fact that deterministic initial 
conditions can be considered more effectively by the POD. The statistical initial conditions 
were drawn from a three, see Eg. (132!) or two, see Eg. (1331) . dimensional subspace which is 
reproduced poorly by a reduction to a 4 dimensional space, which has to consider the time 
evolution also. 
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C. Nonlinear Diffusion 



This system is defined by Eq. (fTT|) . As initial conditions we have chosen a front with 
uniformly distributed position X and normally distributed height H: 

^t = 0,Xi) = —tanh{{xi- X)10). (33) 

Under this conditions all methods were able to reproduce the dynamics well, see Fig. [91 
Surprisingly the full POD method gave poorer results than even the Fourier mode based 
reduction. This is to a lower extent also true for the initialization run of the DMRG POD. 
The iteration lead to an improvement although the 2nd iteration gave similar results as the 
initialization. Further iterations again increase the accuracy, so no general statement can be 
made. After a fast saturation by applying the iteration procedure we observed repeatedly 
a decay in the quality of the results, which we attributeto the accumulation of numerical 
errors. 



V. DISCUSSION 



A. Computational Load 



For all calculation steps, e.g. diagonalization, Gram-Schmidt orthonormalization etc., 

nn 

standard algorithms were applied [6|, |21|]. The focus was more on a concise assessment of the 
new algorithm instead of an optimal solution of the toy problems. For the diagonalization of 



the covariance matrix, e.g. first a Householder-tri-diagonalization was performed [2l|, which 
is an 0{N^) algorithm. The calculation of the POD, either for the complete system or for 
the superblock system was performed with the same routine. This comprised the simulation 
as well as the diagonalization. 

For a POD the simulation of the system in the time-span of interest is additionally neces- 
sary. The required computational load for this simulation is implementation dependent and 
is denoted with S{-). Within our approach the simulation and diagonalization is performed 
only on the superblock system. Comparing the results from Fig. [6] and [7] suggests, that the 
necessary number of iterations (sweeps) does not depend on the full system size A^. If we 
denote the superblock size with M and the number of iterations with Aj a naive estimation 
of the computational load is given in table IV A[ 
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full system POD 


DMRG POD 


0(iV3) + S{N) 


NNiO{M) + S{M) 



TABLE I: Naive estimation of the computational load with full system size A^, superblock size M 
and numbers of iterations Aj. The computational load for the simulation is denoted with S{-) 

For a more quantitative analysis we have measured the time necessary to perform a full 
POD comprised of simulation and diagonalization. Then we did the same for the initializa- 
tion of the DMRG POD algorithm including all simulation and diagonalization steps until 
the superblock system described the full system of dimensionality A^, compare Fig. [21 and 
a first reduced basis had been calculated. We also measured the computing time for one 
further iteration step in the same way as for the initialization. The computing time is con- 
stant for all iteration steps so further data was extrapolated. The underlying equation was 
the deterministic initialized Burgers equation although the choice for an equation affects the 
computational load only marginally. As parameters we have chosen ht = 0.005, d = 0.01 and 
u = 0.1. Fig. [To] shows a logarithmic plot of the results. The DMRG POD approach shows a 
lower amount of computer time for the initialization step. For higher system size this holds 
also for the iterations. Generally the scaling with A^ is favorable. Note, that here only the 
DMRG method should be assessed. For this purpose public assessable standard algorithms 
are sufficient, although much more effective methods could be possible. All calculations were 
performed on an Intel Dual Core machine, using a single CPU. 



B. Stability 

Many numerical schemes and the explicit Euler method in particular show instabilities 
for certain parameter ranges. For the explicit Euler method the stability condition is 

\1 + Xht\<l, (34) 

where A is the largest eigenvalue of the generator of evolution and ht the size of the time step. 
For the Laplace operator the highest frequency component is the first to become instable 
while increasing ht- The eigenvalue is ^ sin^ (^^^^^^^) ~ (^^- Consequently, we should 
have ht < -^^^j d being the diffusion constant. We have performed calculation directly 
at this limit, see Fig. [TT] and have seen no signs of instability. Although the nonlinear 
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FIG. 11: Reduced Burgers equation L^-Error E{t), statistical initial condition, N=40, ht = 0.00625, 
d = 0.05 and u = 0.1. The error is expressed in units of <I>, for the time axis arbitrary units are 
employed. Note that for clarity not all data points are shown as symbols. 
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FIG. 12: Pictorial representation of the approximation. 



Burgers equation was considered, the previous point holds, since the symmetric derivative 
operator has no nonzero real part eigenvalues. By increasing ht the instability appears for 
all methods. These calculations were only first tests and further work is required. However, 
it can be assumed that in the DMRG POD calculation the instability originates only from 
the newly inserted nodes which correspond to the highest spatial resolution. Combined 



implicit-explicit methods 



22j could be used to solve this problem. 



C. Interpretation of the Algorithm 

The various steps, necessary for a DMRG version of the Proper Orthogonal Decomposi- 
tion, seem to be complex on the first glance. It is also not clear why this approach should 
be effective. We now shortly depict the basic idea behind this algorithm. 

The DMRG algorithm decomposes the spatial domain but considers an interaction of the 
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domains due to the superblock concept. The single blocks, together with information about 
interaction with neighboring blocks and the reduced operators, describe spatial regions with 
a higher number of nodes than the number of dofs actually retained in the block. Inserting 
nodes described by the full (yet already discretized) dynamics corresponds to increasing 
the spatial resolution locally. The surrounding blocks simulate the environment for a small 
subsystem with correct dynamics. In the DMRG POD algorithm now the area with high 
resolution is moved through the system. Thereby the parameters of superblock system, i.e. 
the block basis and operator matrix elements, are adapted to approximate the full system. 
For more graphical illustration see Fig. [121 These arguments are a bit heuristic, but until 
now no rigorous proof for the algorithm has be given. 



D. Conclusion and Outlook 



To summarize, we have given a demonstration of applicability for a new algorithm to 
calculate an approximate POD without ever simulating the full system. Our approach also 
makes practically no assumption on the equations that define the dynamics. The approach 
has been tested for linear systems where its performance was even higher than the full 
system POD results but considerable worse than the optimal reduction. Several nonlinear 
systems have been considered. For the Burgers equation the results of the full POD and our 
algorithm were comparable and significantly better than a Fourier mode based reduction. 

Further work on this topic will include extensions to higher dimensional systems. A 
method for 2 and 3 dimensions is currently in progress. Further, driven systems and systems 
with noise shall be implemented. A closer analysis of the quality of the approximations as 
well as the limitations of the method has to be performed. We have access to the amount 
of discarded information from the discarded eigenvalues in the truncations as well as in the 
POD steps. This suggests an adaptive approach for the reduction. 

I would like to thank Prof.F. Schmid, Prof. Ph.Blanchard and Javier Rodriguez- 
Laguna 23|| for discussion and the German science foundation (DFG) for support (DE 
896/l-(l,2) 2004-2007). 
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VI. APPENDIX 



For completeness, we assess in the following the error of the reduced evolution for the 
linear case. For the optimal reduction we require a minimal L^-error for the reduced field 
with respect to the un-reduced evolution. The full time evolution in the dimensional 
phase space is generated by L as 

<|.(t) = e(*-*o)^$(to). (35) 

The explicit Euler algorithm approximates this by 

^t)^{l + htL)^t^). (36) 

We assume that all eigenvalues of L are negative or zero. A positive eigenvalue would lead to 
an unbounded exponential growth in Eq. fl35p which is unphysical. Considering only linear 
projections the reduction is defined by the operator P which is the orthogonal projection 
to the relevant subspace Range(P). P can be constructed from an ONB of this space. 
Equivalently, it can be defined via the ONB (namely C) of Kern(P) so that P = 1 — CC^. 
The reduced time evolution becomes 

<|(t) = e(*-*'')^^^P<l>(to) = e(*-*«)^<l(to), (37) 

since after each (infinitesimal) time step the component within the irrelevant subspace, i.e. 
Kern(P) are projected out. For a general P the eigenvectors of L are not the same as for 
L, but known eigenvectors of L are always the column vectors of C. 



A. Long Time Optimised Projection 

If we assume that the eigenvalues of L are < 0, for long times t ^ 1 the time evolution 
operators e*^,e*^ become the projectors onto the kernels of L or L, respectively. In the 
eigenbasis ipf"^ it is simply 

^f^e^^^f ^ = %e*^' (38) 

The product of the reduced evolution operator e*^^^ and P converges for long times to the 
projector onto Kern(PLP) fl Range(P). More explicitly this is 

^hme*^ = l|Kern(L), (39) 
hme*^^^ = llKern(PLP)- (40) 
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This gives for the error 

= E^^'^(^) = (l|Kern(L) ' l|Kern(PLP) P) "^(^o)- (41) 
In the long time hmit we can obtain a zero error for all initial conditions if we have 

Kern(PLP) n Range(P) = Kern(L). (42) 
This is achieved by requiring 

Kern(L) C Range(P), and (43) 
Range(P) L invariant. (44) 

as we show now. 

Consider a G Range(P). Then Pcj) = (f) and due to the L-invariance of Range(P) it is 
L(f) G Range(P) resulting in PLPcj) = PLcj) = L(p. This gives for P with Range(P) being 
L- invariant 

Kern(PLP) n Range(P) = Kern(L) n Range(P). (45) 

Eq. (l42l) can be retrieved from Eq.( l45l) just by requiring condition ( l43l) . Thus, in the long 
time limit Eq. fl41l) becomes identically zero. 



B. Short Time Optimised Projection 

For short times we consider here the reduction from a A^- dimensional to a (A^ ~ 1)- 
dimensional system. For further reductions the results can be applied by iteration. The 
projector P becomes then Pij = lij — CiCj where c is the removed state. In order to minimize 
the error for the short time evolution measured by the L^-norm we have to minimize 

E,{t) = ||e*^0-e^^^P</>||2 (46) 

^ ||(l+tL-P-PLP)(/)||2 = ||P0||2. 

Here we have already used an expansion in powers of t and truncated after the first order 
terms. 

Since we have no information on (j), we minimize Eq.( H6l) by using the Frobenius norm 



2l|, i.e. 



\f of the error operator E. The Frobenius norm is consistent with the L^-norm 

ll^xlla < |A|f||x||2 Vy4 G P"^",x G P". (47) 
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By inserting F = 1 — C we get for the error operator 



l-P + t{L-{t- C)L{t - C)) 
= C + t{L - L + LC + CL - CLC) 

C + t{LC + CL-CLC). (48) 

We assume L to be symmetric, i.e. Lij — Lji. Thus L has an orthonormal eigenbasis 

{^ia}a=i N where the columns arc the eigenvectors of L. The eigenvalues are and the 
matrix elements of the error operator E are decomposed in this basis as 



Cij ~\~ t ^ ^ j LinCnj ~l~ CijiLyij ^ ^ CijiLjimC'mj j 
n \ m / 



(49) 



with 



Qj = '^<^aiCaCi3<^l3j, (50) 
a/3 

mn aPnm 

LjnCnj = 'PanLinCgCisippj , (52) 

(PaiCaCpLnjippn- (53) 

n a/3n 

We use the orthogonality of the <^q,, i-e. 

^2 Vo^i^Pi = = ^ <fia'Pil3- (54) 
i j 

In the eigenbasis the removed degree of freedom c can be written as c with components 

Cj = ^ (PaiCa , = ^ Vai^PiCa = ^ VfiiCi- (55) 

q; ia i 

The average of L in the removed state c, is 

(-^)c ■ ^ ^ ^nLnmCm ^ ^ 'finiC-iLnm^mjCj (56) 

nm nmij 
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The matrix elements from Eq.s fl50H53ll become 

Cij = CjCj, (^'^) 

^ ] CjnLnmCmj = Cj(L)^Cj, (58) 
mn 

n an 
n 

= (59) 

n n 

= CiCjXj. (60) 
Thus for the matrix elements of the error operator we obtain 

E,, = c~c, (1 + t (A, + A, - (L) J) . (61) 
We minimize the Frobenius norm of E given by 

\e\f = y. = E s'^l (1 + ^ (^^ + - Wc))' (62) 

for a normalized c, i.e. 

l = \\c\\l = J2^' = J2^'- (63) 

i i 

Since is a linear operator it follows that ||i?x||2 = | |x| ^ |-Ex| I2 with x = x||x||2. Without 
no restriction to ||x||2 the zero vector would always minimize ||i?x||2. Furthermore, each 

lead to the same x with llxlU = K. This is not true for the 



lower bound K for ||x||2 wil 
general nonlinear case as in 



24|. 



Incorporating this condition \E\p reduces to 

\E\l = l + 2t (L), + 2t (L), - 2t (L), + {L)l 

+ 2^ {L)l - 2t^ {L)l + e {L)l - 2^ {L)l + f {L)l 

= l + 2t{L)^ + t'{L)l = {l + t{L)S 

^ \E\f = \l + t{L)J. (64) 

Consequently, in order to minimize \E\f we have to minimize (L)^. 



25 



The minimization itself is performed using Lagrangian multipliers for the constraint 
Eq. fl63l) . The necessary condition for a minimum is 

= J: E (5?^' + -' (!-«?)) (65) 

i 

= 2ck (Afe - 1]) . 

This is true if either = or = A^. The last equation can only be true for a single value 
of Afc. We denote the nonzero component as Ck' 7^ and = 5kk'Ck'- From Eq. fl63|) it follows 
further that Ck = Skk'- 

Inserting this in Eq. llMl) we obtain 



\E\ 



k 



\l + tXk'\. (66) 



For small t , i.e. t < |Aj|~^ this is clearly minimal if we choose A^' to be the smallest 
eigenvalue. 

Further iterations, e.g. n times, of selecting the irrelevant states remove successively the 
eigenstates corresponding to the n lowest eigenvalues. This is due to the fact that the spaces 
Kern(C) = Range(P) and Range(C) = Kern(P) are by construction L invariant. This also 
makes the iteration unambiguous, a feature that is in general not present for nonlinear 
problems. 

Note also that since Aj < the reduced states always belong to Range (L) as long as any 
remaining eigenvalue, i.e. an eigenvalue of P„_iLP„_i is nonzero. Here, P„_i results from 
the previous reduction step. In this case the error always vanishes for long times. 

Summarizing, the optimal short time projection leads to results that are not only consis- 
tent with the long time accuracy requirements, but even include them. 
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